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p^j \ In many dynamical systems there is a large separation of time scales between typical 

events and "rare" events which can be the cases of interest. Rare-event rates are quite 
difficult to compute numerically, but they are of considerable practical importance 






e : in many fields: for example transition times in chemical physics and extinction times 

c^ ■ in epidemiology can be very long, but are quite important. We present a very fast 

numerical technique that can be used to find long transition times (very small rates) 

Cd ■ 

^ ■ in low-dimensional systems, even if they lack detailed balance. We illustrate the 

'-^ . method for a bistable non-equilibrium system introduced by Maier and Stein and a 

a ■ 

O ' two-dimensional (in parameter space) epidemiology model. 
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I. INTRODUCTION 

Important events for a transition may have a time scale many orders of magnitude larger 
than typical events; in this case, they are called "rare events." They have been studied 
in a number of different contexts, including the extinction of diseasea^ or of populations^, 
network queue overflow^, and slow chemical reactions^. 

To fix our ideas we consider a very simple problem, that of an endemic disease which 
fluctuates to extinction. Consider a population of fixed size, N, with S members who are 
susceptible to an infection, and I who are infected. When S encounters /, the infection is 
transferred with rate PSI/N; /3 measures the infectivity. Infected agents can spontaneously 
recover with rate k, and can be immediately reinfected. This is called the SIS (Susceptible- 
Infected-Susceptible) model- . In this simple form it can be thought of as a Markov process 
in the number of infected; note that S = N — I. Thus: 

W(I -^ / + 1) = PSI/N 

(1) 
W{I ^ I-1) = kI. 

There is a long-lived state where the disease persists when Rq = P/k, > 1, namely I = 
N{1 — Rq^)- However, there is an important rare event, namely a fluctuation to / = 0, 
which means that the disease is extinct and cannot return. For this simple form of the 
model the mean exit time for this transition, T, can be found exactly^. For large A^, there 
is an asymptotic formula: 

r -V ^0 1^^ NdofrRo-l+l/Ro) (r,\ 

' ^(/?o-l)^Viv" • ^'' 

This formula has features that are generic to the kind of "barrier climbing" problems that we 
treat here. The exit time is of the form g{N) exp{NW) where g is a slowly varying prefactor, 
and ly is a generalized barrier height (or quasipotential) scaled by the large parameter, A^. A 
plot of the exact results for the SIS model with _Ro = 2 is given in Figure [1] along with results 
from numerical computations that we will describe below. Note that even for modest-sized 
systems the mean time to extinction can be huge: for A^ = 300 we have T ^ 4 x 10^^. 

As we see from this simple case, these phenomena are frequently out of reach for brute- 
force simulations. To overcome this problem, many techniques have been developed^. In 
this article, we revisit this problem and present a very efficient technique which we call the 
barrier method and which gives the mean first-passage time for a transition to an unlikely 
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FIG. 1. The mean time to extinction, T, as a function of the population size, N, for the simple 
SIS model calculated exactly^ and by two of the numerical methods to be described in the text. 
All numerical results agree with the exact result within statistical error. 

target state. The method does not depend on special features such as knowledge of a steady- 
state distribution or detailed balance in the process. We need only that the dynamics be 
stochastic, and that the transition probabilities depend on the current state. That is, we 
deal with Markov dynamics, which can be reversible or irreversible. 

The essence of the algorithm is to follow the development of an ensemble of systems and 
oversample the cases that happen to approach the target, and not allow backtracking away 
from the target. In this respect our method resembles the signposting algorithm^i^ that we 
developed for finding the penetration of a random walker into a fractal. 

We believe that the barrier method is the most efficient method available for computations 
in relatively low dimensions. In this paper we explain how it differs from previous techniques, 
and we apply it to two systems. The first is described by a non-equilibrium model introduced 
by Maier and Stein.— The second is a generalization the SIS model in which the population 
is allowed to fluctuate^. In this model, we are concerned with the average time for the 
disease to go extinct. Lastly, we discuss the advantages of our method and future work. 

II. BACKGROUND 



To study rare events in a Markov process we must deal with states in state-space that are 
unlikely to be visited in any simulation of reasonable lengtb^^. Most sample paths spend the 



majority of time visiting the most likely states and give good estimates of the corresponding 
probability density. For rarely visited regions, we need to use special methods. 

An example of such a method is biased sampling, in which is we arrange our simulation 
to be biased towards rarely visited regions. Two important subsets of this approach are im- 
portance sampling and splitting techniques. Importance sampling, the most commonly used 
method for equilibrium systems, requires some a priori information about the probability 
distribution. In contrast, if the probability distribution is only accessible via simulation, 
splitting and related techniques are very useful. The latter case is the focus of our work. 

A. Splitting and RESTART 

Splittingi^ involves placing a barrier in state space. When a sample path crosses the bar- 
rier it is split into independent realizations whose statistical weights add up to the original. 
By placing one of these splitting barriers in a region which would be infrequently visited, 
that region will subsequently be better sampled. For very difficult-to-reach regions of state 
space one barrier is not sufficient. The use of several barriers is called multilevel splittingi^. 

Multilevel splitting has two drawbacks. First, if the barriers are too close together or too 
far apart, the number of simulations will grow or decay exponentially. Further, realizations 
which have small weight (because they have been split many times) will often backtrack, 
i.e., tend to move back to the well-sampled regions and waste computational effort. Some 
of these problems have been solved by RESTART= (REpetitive Simulation Trials After 
Reaching Threshold) which is designed for dealing with queuing problems. 

In RESTART^* one considers nested subsets of phase space, A D B D C D D, where 
D is the target. Barriers are placed between these regions. A sample path is started in A 
and evolves until it reaches B. Then the sample is split into R 'retrials' with equal weight. 
One of these is designated the primary, and all realizations evolve independently. If any of 
the non-primary paths backtracks into A, it is terminated. Each barrier is crossed in turn, 
and the time spent in D by the reweighted paths gives an estimator of the phase space 
probability in D. 

RESTART partially solves the backtracking problem because most samples that exit low 
probability regions are terminated. However, in the original version^^ii^ it still can lead to 
a divergence in the number of samples if the barriers are too closely spaced. It has been 



noted that the barrier placement problem could be partially alleviated by performing fixed 
effort RESTART instead of fixed splitting RESTART^, but this does not completely fix the 
problem. We give another approach to this problem below. 

B. Forward Flux Sampling 

Forward fiux sampling (FFS)-^ uses the same principles as splitting in applications to 
computational chemical physics. Most rare-event techniques in this area^ require equilibrium 
ensembles and detailed balance. FFS is unusual because it is applicable to systems without 
detailed balance. It has been used to study genetic switche3i2,-i£^ nucleationi2,"Ji and a model 
problem due to Maier and Steini^i^ which we also treat below. 

FFS finds the first passage time between metastable states A and D as follows. First, we 
run a single long simulation and count the number of times the sample path exits A through 
barrier Aq, which bounds region A. The average fiux through the barrier, ko, is gotten 
by dividing this number by the total simulation time, discounting the time associated with 
trajectories that reach D and return to A. If we call Xm the barrier around D, the transition 
rate from A to D is: 

kAD = koP{XM\Xo)- (3) 

where P(Am|Ao) is the probability that a sample path which starts on Aq will cross Am before 
going back to A. 

We can get P(Aa/|Ao) efficiently by introducing intermediate barriers Xi,i = 1 . . .M — 
1 divide the sample space along level surfaces of some reasonable guess for the reaction 
coordinate - we call this the order parameter. The probability factors into: 

Af-l 

P{Xm\Xo) = n Pi^^+l\^^)^ (4) 

i=0 

where P(Aj+i|Aj) is the probability of starting at Aj and reaching Aj+i before going back to 
Aq. To measure P(Ai|Ao), R samples are started from the locations along Aq where they left 
A in the first step. The paths are continued until they reach Ai or go back inside Aq. The 
fraction of that reach Ai is the estimator of P(Ai|Ao). Then we proceed to A2 and start R 
samples, etc. The point is to break down a long sample path into a series of short segments. 
FFS does not allow the number of samples to diverge, as in splitting. However, it does 
backtrack because samples which start at Aj must be allowed to return to A. This effect 



can be somewhat reduced by pruning of the backtracking paths; see^^. However, if there are 
metastable states the region between A and D, backtracking can take a long time. 

Further, the calculation of k^ requires that the initial long simulation reaches the end 
state D at least once in order to properly sample the entire region between the start and 
the target. If one does this, FFS will often be bottlenecked by the calculation of ko; this 
defeats the purpose of using a rare-event technique. Fortunately, for systems with featureless 
barriers, running an initial simulation that crosses Aq a fixed number of times, say lOR, is 
typically sufficient. We call this version "approximate FFS." For a comparison of the two 
versions of FFS for the SIS model, see Figures [H El The approximate method of FFS gives 
good accuracy for this problem, and is quite fast. However, for systems with metastable states 
in the region between A and D, this method of determining ko is not sufficient. Our method 
(see below) overcomes all of these problems. 

C. Milestoning 

Milestoning^^ is a technique for equilibrium systems in which one runs simulations of 
short paths between barriers to find the local first passage time from one barrier to the 
next. The equilibrium ensemble on the barriers gives the launching points, and there is no 
backtracking at all. The local first passage times are put into an integral equation to find 
the global first passage time. As we will see, we avoid using equilibrium considerations on 
the barriers by keeping track of the landing points of individual paths, but otherwise, our 
method uses similar ideas. 

III. BARRIER METHOD 

A. Algorithm 

In the barrier method, we consider the same sort of problem as for FFS. For the moment, 
we assume that the locations of the barriers are known a priori, as shown in Fig. |2l We will 
discuss the best way to distribute the barriers below. In the first step, R trial simulations 
are started in A and run until they reach Aq. Each trial r ends at W^^ at time r^ )^^^. R more 
trials are started at the locations along Aq where each r ended in the previous step and run 
until they reach Ai. These trials can backtrack as far as they need to. The locations along 




FIG. 2. Barrier method. Three different paths are started from A. The paths cross barrier then 
cross barrier 1. The middle path then moves backwards and reaches barrier 0. The end points on 
barriers and 1 are used to "jump" the path back to the barrier 1, grey path. Paths terminate 
once they reach B. 



Ai where the trials stopped are W{_^ and the transition times from Aq to Ai are t^^ ^_^. For 
each r the total time to start from A and reach Ai is r^ ^^^ = ''"a A) + '^x Xi- 

On the next step we eliminate backtracking. We start each trial r at the location on Ai 
where the previous trial r finished, W^ . Each sample path continues until it either reaches 
A2 and stops or returns to Aq. If the trial returns to Aq, we move it back to Ai and add a 
sample of the tim,e to return from Aq to Ai. In practice, we find the the closest W^^^ to the 
current location and add r^^^ ^^ to the trial time and continue the trial at W^^ . That is, we 
move the sample point "in one step" to the next barrier. Continuing, the sample path can 
either reach Ai and stop or go back to Aq, in which case we repeat the process of jumping. 
The locations along A2 and the transition time for each r are VFJ^ ^^"^ ''"Ii X2- The total time 
to reach A2 is t\ ^^ = ''"a Ai + '^Xi X2 ■ Then we repeat the process for the next barrier, and 
continue until r^ ^ and Wam ^^^ calculated. The average transition time from A to Xm is 
^r-rkxjR- 

This method differs from FFS in two ways. First we work with transition times instead 
of transition rates and never deal with probabilities directly. Second, the barrier method 



does not require sample paths to travel from Aj all the way back to Aq. This can lead to a 
dramatic improvement in efficiency over FFS, as we will see in the next section. Note that 
metastable states in the barrier region pose no problem for this method. 

B. Accuracy and Efficiency 

The simple SIS model described above is exactly soluble. This allows us to make a direct 
comparison between FFS and the barrier method. Using A^ = 20 to 200, we found less 
than 1% difference between the barrier method and the exact results, which was within the 
variance of the measurement. Both versions of FFS gave similar accuracy; see Fig. [1] The 
exact version of FFS, which samples the whole region between A and D in the estimate of 
ko is impractical for A^ > 60, as shown in Fig. [31 

We now compare the efficiency of approximate FFS and the barrier method. Simulation 
efficiency can be defined by the amount of computation time, C, needed to obtain the exit 
time within a relative accuracy of a; we chose a = 0.1. Using this definition, we varied 
the number of trials, R, to find C for Rq = 2 for a range of A^, from 50 to 290 in steps of 
10 for FFS and the barrier method. For both we used A^/10 — 1 barriers, placed every 5 
infectious population size steps starting at / = 5. No attempt was made to optimize the 
barrier placement. The comparison is given in Figure [31 The simulations were performed 
on a 2.8 GHz Intel processor. We find that both algorithms appear to have a power-law 
relationship between C and A^, with powers of roughly 3 and 2 for FFS and the barrier 
method respectively. This shows that for a case where the barrier method is not obviously 
better, because of the simple landscape, it still significantly outperforms FFS. 

C. Dynamic Barrier Placement 

We have also developed a way to determine the locations of all of the Aj barriers as the 
simulation progresses. The method outlined in this section can be applied to FFS as well as 
any other biased sampling methods that use barriers, e.g. RESTART^^. We have previously 
used a similar method for signposting^i^. 

Assuming that the location of the ffist barrier has been chosen and all of the W^ values 
have been calculated, p probe trials are started along Aq at W^^ for some randomly chosen 
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FIG. 3. a.) CPU time required to calculate T to a fixed precision (arbitrary units) for the SIS 
model. For FFS and the barrier method, the same number of barriers were used and placed in 
the same locations (every 5 steps in number of infected), b) Efficiency comparison between the 
(approximate) FFS method and the barrier method. The solid lines represent the best fit to a 
power-law: with powers 3.3 (FFS) and 2.1 (barrier method). These values should be taken as 
rough estimates of the powers. 

r's. The probe trials are run for some fixed number of simulation steps, S. The location 
of Ai is chosen as at the average longest excursion along the order parameter of the probe 
trials. Next, W^^ and rj^^^ y^^ are recorded. We repeat the process by running more probe 
trials for 5* simulation steps, not counting any steps involved with jumping from Aq to Ai, 
starting at random W^.^s. The average furthest excursion is used as the location of A2. This 
process is repeated until Am is reached. An alternative approach would be estimate the 
location of the next barrier such that a fixed fraction of samples reach Aj+i before reaching 
Aj_i when starting at Aj. An rough estimate of the optimal forward fraction is^^ e~^. We 
did not chose this approach because we found that for small noise the probability of making 
any forward progress is significantly smaller than e~^. This would cause the barrier method 
to stop making progress. 

There has been some previous work on barrier placement or staging using FFS^. The 
scheme uses two FFS calculations: one with a guess of the best barrier locations and fixed 
i?, and a second with either optimized barrier locations or optimized the number of trials R 
for each barrier. The authors found that it is better to optimize the spacing of the barrier 




FIG. 4. An example of the barrier method when the most probable exit path (MPEP) crosses 
barriers several times. The line from A to B represents the MPEP. The x's represent the initial 
locations where simulations reach a given barrier. These locations are used to generate the original 
set of lookups. The diamonds represent the new lookups generated by trials originated on the 
second barrier. Similarly, the triangles are new lookups generated by the trials started on the third 
barrier. Note that without adding new lookups, the simulation would never finish because the 
closest lookups on previous barriers would move the simulations backwards along the MPEP when 
they are moved to the next barrier. 

to get uniform p(Aj+i|Ai) with fixed R than to optimize R for each barrier. When we apply 
dynamic barrier placement to FFS, we also take this view. The important difference from 
the approach outlined above is that we determine the location of the barriers one at a time 
instead of all at once. This leads to computational gains when the trial barrier placement 
is significantly different from the optimal placement. 

D. Additional lookups 



As with FFS, the barriers are should be set up at fixed values of the order parameter. 
The closer the order parameter is to the true reaction coordinate, the more efficient the 
method becomes. For the barrier method, an additional issue arises if the most probable 
exit path (MPEP) crosses a barrier more than once, so that it goes backwards. The original 
set of lookups move simulations from a previous barrier to the current barrier along the paths 
already discovered. These lookups cannot move the simulations along the correct MPEP, as 
seen in Fig. HJ This problem can be fixed by allowing additional lookups to be created if 
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needed. For the barrier method apphed to the generahzed SIS model, below, we used this 
feature. 

IV. MAIER-STEIN MODEL 

Maier and Stein^^ introduced an interesting example of a dynamic system which lacks 
detailed balance. It has received considerable theoretical^^i^^, experimental^i'^^, and compu- 
tational interes1r^'2^'2Ii2£. We study two aspects of this model: the mean exit time, T, from 
one of the metastable states and the distribution of exit locations along the separatrix. 

The model is specified by two coupled stochastic differential equations: 

x = /^(x)+^^(t), 

(5) 

where x = {x,y) and f = {fx, fy) is the time-independent drift field: 

fx=x — x^ — axi? , 

(6) 

fy = -fiyil + x"^). 
For a = fj, the model obeys detailed balance. The white noise ^ = (^x,Cy) has variance e: 

im) = 0, m + r)m) = ^^^a^ - ^)- (7) 

We are interested in the small noise case: e — )■ 0. The model is bistable with the metastable 
states located at x = (±1,0). There is a separatrix at x = 0. The exit time, the transition 
time from one of the metastable states to the separatrix, is, in the small noise limit: 

T oc e^/^ (8) 

where W is the generalized barrier height. In this model, W is a function of a and /i. If 
H = 1, for 1 < a < 4, there is a unique MPEP, and for a > 4 there are two MPEPs. 

A. Simulation 

To simulate the system, we solve Eqs. (|5]) using the Euler method: 

x{t + h) = x{t) + h [x{t) - x{tf - ax{t)y{tf) + \/lh, 
yit + h) = y{t) + h {-fiy{t){l + ^(t)^)) + Vdl, 
where h is the time step. 
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FIG. 5. Barrier height, W, for the Maier-Stein model as a function of a for // = 1. Circles, barrier 
method; line, analytic theoryi^ for e — )> 0. 

B. Exit time and distribution of exit points 

We start R simulations at the left metastable state (—1, 0) and iterate Eq. (^ to obtain 
the trajectories. We dynamically locate each barrier. The first barrier is located at the aver- 
age furthest excursion along x, the order parameter, using 5/4 probe steps. All subsequent 
barriers are located after S steps. To calculate the exit time, we have three parameters 
at our control: the number of trials, R, the time step, h, and the spacing of the barriers, 
controlled by S. To find W, we record the exit time for various values of e and find the slope 
of of 1/e vs. In T. 

C. Results 



Using the barrier method we measured W and the exit distributions along the separatrix, 
P{y). We computed W^ for yU = 1 for a = 1 to a = 8, and compared to analytic theoryi^. 
This results are shown in Fig. O 

In order to use this model, we have to choose h, R, and S. These parameters have 
different, competing effects. Increasing h increases the efficiency but decreases the accuracy. 
The opposite is true for S. Also, effects of the values of h and S are connected; small S 
causes there to be a large number of barriers, which requires a small h to give an accurate 
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FIG. 6. The exit distribution for the Maier-Stein model for a = 10, ^ = 0.67 and e = 
0.0025,0.0004,0.0001. The simulation results are denoted by symbols. The dashed line is the 
symmetrized Weibull distribution P{y) = A^lyp'^"^ exp (— |y/Ap^/e) from theor y^^'^*^ , where N 
is the normalization, A is a parameter of order unity. Solid line: Weibull distribution convolved 
with a Gaussian with a = Byfe. 

result, and vice versa. We used h = 10^^, R = 7 ■ 10^, and 5 = 4- 10^. However, we used a 
value of S which is four times smaller than the general case to locate the first barrier. The 
exit time was measured with five independent trials for five values of 1/e: 20, 40, 60, 80, 
and 100, for each a. The calculated value of W and the small-noise theory^^ are shown in 
Fig. [5l The numerical results are consistently smaller than the theory and increase as 1/e 
increases. The values of W in Fig. [5] were consistent over other values of h, R, and S, that 
we tested. 

We ran separate simulations to measure the exit distribution along the separatrix. Be- 
cause we are interested only in the final location along the separatrix and not the time it 
takes to reach it, we found that we could use significantly larger h. We were able to get reli- 
able results for values of h as large as 0.0001. This allowed us to reach much smaller values 
of e; see Fig. El The parameter values used to obtain the results in Fig. |6]are: h = 0.0001, 
R = 10^, and 5* = 4 ■ 10^. The results were averaged over twenty independent simulations 
for each e. The significant result is that for the smallest noise value, e = 0.0001, the value of 
P(0) is not close to zero as the theory suggests. Rather the ratio of -P(O) to the maximum 
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of P{y) appears instead to be increasing as e decreases. 

These results are consistent with other simulations^ and experiment al^S results. Previ- 
ously, the results were assumed not to match the theory because the values of e were not 
small enough. The barrier method allows us to reach a value of e which is 50 times smaller 
than the best previous simulation result^^ and 110 times smaller than the best experimental 
result^ without P{0)/Pmax getting any closer to zero. 

We propose that the reason for this is that the theoretical prediction of the Weibull dis- 
tribution represents the leading term in e. For finite e the distribution should be "rounded" 
over a scale y = 0{^/eY^. Accordingly, we convolved the asymptotic theory, the Weibull 
distribution given in the caption of Fig. [6l with a Gaussian with a = B^\ see Fig. [6l We 
find good agreement between the simulation results and the convolved theory. The values 
of B which gave the best fit were 0.8, 0.85, and 0.85 for e = 0.0025, 0.0004, and 0.0001, 
respectively. The fact that B is roughly constant over a wide range of e gives support to our 
estimate. 

The reason that P(0)/Pmaa; does not tend to zero is that even though the rounding is 
over a scale that decreases as -\/e, the location of the maxima of P{y) also moves toward the 
origin. These locations are y^ax = ±2"''/^Ae'^/^(2 — ;u)^/^, so that \ymax\ «: Ae^^^"^ ^ Ae^^^ 
for n = 0.67. If A were constant, P{0)/P{ymax) should approach zero rather slowly. We 
numerically found this rate to be e°'^. However, we find that A is not constant. Our 
best fit values for A, which is unaffected by the convolution, are 0.94, 0.72, and 0.4 for 
e = 0.0025, 0.0004, and 0.0001, respectively, so that P{0)/P{ymax) does not approach in 
our computations. 



V. GENERALIZED SIS MODEL 

It is interesting to generalize the SIS model of ([T]) to allow fluctuations of the total 
population by introducing birth and death rates^'i^'^'^. Now there are two independent 
stochastic variables, S, the number of susceptibles, and J, the number of infected, and four 
parameters: /i, the birth and death rate assumed equal, /3, the infectious contact rate, k, 
the disease-recovery rate, and A^, the steady-state population size. The transition rates are 
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now: 

W[iS,I)^iS + l,I)]=fiN, W[iS,I)^iS-l,I)] = fiS, W[iS,I)-^iS,I-l)]=fiI, 
W[{S,I)-^{S+1,I-1)] = kI, W[{S,I)^{S-l,I + l)]=f3SI/N. (10) 

The model has an endemic state when i?o = /3/(/i + k) > 1. There is one stable fixed point, 
the endemic state {S, I) = (NRq^, N(1 — Rq^)), and an unstable saddle point (A^, 0), where 
the disease is extinct. We seek the transition time from the endemic state to the disease-free 
state which will be of the form T ~ exp{NW) , as above. 

We will be interested in the case of small /i so that population fiuctuations are slow 
compared to disease dynamics. It might seem that the situation would be very similar to 
the case /i = treated above. However, this is not true^^. Population fiuctuations make 
extinction of the disease much easier: the most likely exit path is via a population decrease 
at fixed S followed by extinction along a path of smaller fixed population, and then an 
increase of population of susceptible individuals to S' = A^. 

A. Simulation 

We simulate the SIS system using standard techniques^^'^. We define the order parameter 
as —J; the barriers are added at decreasing values of /. Because this system is on a 2d 
parameter lattice it is much easier to dynamically add new lookups if previously unexplored 
backward regions are reached or if a given lookup has been used too frequently. 

The method we use to add lookups is as follows. Every time a sample moves backwards 
and reaches the previous barrier, a new lookup for that site is generated with a probability 
Pg = Lg/D{S, I), where Lq is a constant which controls the growth rate of the lookups and 
D{S, I) is the number of lookups at site (5, J); if pc is greater than 1 a new lookup is always 
added. If a new lookup value is needed, the value is produced by starting a path at {S, I) 
until it reaches the next barrier, and stops or reaches the previous barrier. There a lookup 
is used to move back to the current barrier to continue. Note that a new lookup can move 
back to a previous barrier and cause another new lookup to be generated at that previous 
barrier. This cascading effect can continue until the first barrier is reached. This effect 
makes programming the algorithm more complex, but the cascading is necessary to produce 
accurate results when the MPEP crosses some of the barriers several times, as illustrated in 
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FIG. 7. a.) Generalized barrier in the SIS model for /i = 1, k = 100 as a function of Rq. The line is 
the theoretical estimateSS. b.) Distribution of first passage points on barriers for N = 1000, i?o = 
1.667, /i = 0.25. The exit point on the line / = is in good agreement with the predictio'"^ 
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B. Results 

We found the exit time in SIS model, T ~ exp(A^iy), for /i = 1 and n = 100 and varied (3 
to obtain different values of Rq- We compare to analytic work— i^ii^ in Fig. [7] and Fig. |8l In 
these references the authors chose fi = 0.02. They did this because the very large separation 
in time scales allowed them to use their analytic techniques. 

We calculate W by performing a linear regression of In T versus A^, as in the Maier-Stein 
section. We ran 10 simulations for each of following N: 50, 100, 150, 200, 250, 300, 350, 
400, 450, and 500. For the last three values of A^ we did this only for Rq < 2.5. By looking 
at the residual of the fit, we found for smaller Rq the smallest sizes were not large enough 
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FIG. 8. Exit paths for the SIS model for N = 700, Rq = 1.25. Left panels, some sample paths, right 
panel, average over many paths to estimate the MPEP with a dashed line for the corresponding 
theor}*2^. 

to reach a constant value of W. The smallest A^ included in the fits were 350, 350, 300, 300, 
250, 250, 200, 200, and 200 for Rq = 1.5, 1.625, 1.75, 1.875, 2.0, 2.125, 2.25, 2.375, and 2.5 
and greater, respectively. The calculated values of W are shown in Figure [7] and compared 
to analytic estimates^. 

The barrier method does not give the exit path directly. If we plot the distribution of the 
ly 's on the barriers, it is the distribution of first passage points. However, the paths that 
continue are not uniformly distributed on the barriers. Nevertheless, it is interesting to plot 
the first-passage distribution. Fig. [3 

We did a separate, brute force, computation to find the actual MPEP. A few sample 
paths are shown in Figure |8] along with the average of many paths. 

VI. DISCUSSION 

In this paper, we developed a new rare-event technique, the barrier method. We described 
the relationship between it and other related techniques, such as FFS, and showed that for 
a simple model problem the barrier method is more efficient than FFS, especially for large 
N. 

The barrier method was then used to find exit times and distributions for the Maier-Stein 
model^. We found fairly good agreement with theory on exit times and good agreement 
with previous simulation and experimental results, on the exit distributions. We determined 
that convolving the theory with a Gaussian to account for the next correction to the theory 
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gave excellent agreement for the distribution of exit points. 

The exit times for an SIS model with births and deaths were then calculated. The results 
agreed with analytic estimatesiii^i^. The MPEP was also found for this case. 

The barrier method is an excellent tool to determine rare events in low dimensional sys- 
tems. In this paper we have treated one and two degrees of freedom, and we have preliminary 
work for three dimensions. We believe the most important aspect of the barrier method is 
the elimination of practically all 'backtracking.' This can allow traversal of landscapes with 
many metastable states . The method is also general enough to apply to on- and off-lattice 
problems, equilibrium and non-equilibrium problems and any system that can be written as 
a non-deterministic Markov process. 

There are limitations to the method. The most important is that high-dimensional prob- 
lems (i.e., systems with many degrees of freedom) are difficult to treat this way. For example, 
for nucleation problems in Ising systems of A^ spins the number of dimensions is 2^. The 
barrier method eliminates backtracking so that each barrier must be sampled well enough 
for the lookups to be accurate, a daunting task in very high dimensions. Also we need to 
determine which lookup is closest in a high-dimensional space and then figure out if that 
lookup is 'close enough' or whether another lookup must be created. 

However, the point of view that we take, focussing on times rather than rates, is useful 
even in high dimensions. We have also developed a variation of the FFS technique which 
uses this point of view and applied it to a non-equilibrium nucleation problen>2^. 
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